home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 46 / Amiga Format CD46 (1999-10-20)(Future Publishing)(GB)[!][issue 1999-12].iso / -in_the_mag- / reader_requests / scilab / demos / control / lqg.dem < prev    next >
Text File  |  1999-09-16  |  3KB  |  102 lines
  1. exec(SCI+'/demos/control/scheme.dem');
  2. s=poly(0,'s');z=poly(0,'z');
  3. x_message(['Simple example of SISO LQG Design';
  4.            'Demo is in file '+SCI+'/demos/control/lqg.dem';
  5.            'Computes the LQG compensator and plots response'])
  6. n=x_choose(['Continuous time';'Discrete time'],'Select time domain');
  7. select n
  8.  case 0
  9.   return
  10.  case 1
  11.   s=poly(0,'s');
  12.   str='(s-1)/(s^2-5*s+1)';
  13.   rep=x_dialog('Nominal plant?',str)
  14.   if rep==[] then return,end
  15.   Plant=evstr(rep);
  16.   Plant=syslin('c',Plant);
  17.  case 2
  18.   z=poly(0,'z');
  19.   str='(z+1)/(z^2-5*z+2)'
  20.   rep=x_dialog('Nominal plant?',str)
  21.   if rep==[] then return,end
  22.   Plant=evstr(rep);
  23.   Plant=syslin('d',Plant);
  24. end   
  25.  //Nominal Plant
  26. P22=tf2ss(Plant);    //...in state-space form
  27. [ny,nu,nx]=size(P22);
  28. x_message('Now enter weighting matrices');
  29. rep=x_matrix('x-weighting matrix',eye(nx,nx))
  30. if rep==[] then return,end
  31. Qx=evstr(rep);
  32. rep=x_matrix('u-weighting matrix',eye(nu,nu));
  33. if rep==[] then return,end
  34. Qu=evstr(rep);
  35. bigQ=sysdiag(Qx,Qu);
  36. rep=x_matrix('x-noise covariance matrix',eye(nx,nx))
  37. if rep==[] then return,end
  38. Rx=evstr(rep);
  39. rep=x_matrix('y-noise covariance matrix',eye(ny,ny))
  40. if rep==[] then return,end
  41. Ry=evstr(rep);
  42. bigR=sysdiag(Rx,Ry);
  43. [Plqg,r]=lqg2stan(P22,bigQ,bigR);     //LQG pb as a standard problem
  44. Klqg=lqg(Plqg,r);                     //LQG compensator
  45.  
  46. disp(spec(h_cl(Plqg,r,Klqg)),'closed loop eigenvalues:');    //Check internal stability
  47. [Slqg,Rlqg,Tlqg]=sensi(P22,Klqg);  //Sensitivity functions
  48.  
  49. disp(clean(ss2tf(Slqg)),'Sensitivity function');
  50. disp(clean(ss2tf(Tlqg)),'Complementary sensitivity function');
  51.  
  52. x_message('Closed-loop response');
  53. resp=['Frequency response';'Time response'];
  54. while %t do
  55. n=x_choose(resp,'Select response(s)');
  56. select n
  57.  case 0
  58.   disp("LQG demo stops!");break;
  59.  case 1
  60.   xbasc(0);xset("window",0);xselect();bode(Tlqg);
  61.  case 2
  62.   if Plant(4)='c' then
  63.    defv=['0.1','20'];
  64.    title='Enter Sampling period and Tmax';
  65.    rep=x_mdialog(title,['Sampling period?';'Tmax?'],defv)
  66.    if rep==[] then break,end
  67.    dttmax=evstr(rep)
  68.    dt=evstr(dttmax(1));tmax=evstr(dttmax(2));
  69.    t=0:dt/5:tmax;
  70.    n1=x_choose(['Step response?';'Impulse response?'],'Simulation:');
  71.     select n1
  72.       case 1 then
  73.       xbasc(1);xset("window",1);xselect();
  74.       plot2d([t',t'],[(csim('step',t,Tlqg))',ones(t')])
  75.       case 2 then
  76.       xbasc(1);xset("window",1);xselect();
  77.       plot2d([t',t'],[(csim('impul',t,Tlqg))',0*t'])
  78.     end
  79.   end
  80.   if Plant(4)='d' then
  81.    defv=['30'];
  82.    title='Tmax?';
  83.    rep=x_mdialog(title,['Tmax='],defv)
  84.    if rep==[] then break,end
  85.    Tmax=evstr(rep);
  86.    n2=x_choose(['Step response?';'Impulse response?'],'Simulation:');
  87.     select n2
  88.      case 0 then
  89.       break;
  90.      case 1 then
  91.       u=ones(1,Tmax);u(1)=0;
  92.       xbasc(1);xset("window",1);xselect();
  93.       plot2d([(1:Tmax)',(1:Tmax)'],[(dsimul(Tlqg,u))',(ones(1:Tmax)')])
  94.      case 2 then
  95.       u=zeros(1,Tmax);u(1)=1;
  96.       xbasc(1);xset("window",1);xselect();
  97.       plot2d((1:Tmax)',(dsimul(Tlqg,u))')
  98.     end
  99.   end
  100. end
  101. end
  102.